#include "CROP.h"


using namespace std;

void InitCROPParameters(CROPParameters &CParameter)
{
    memset(CParameter.Block,1,sizeof(int)*20);
    CParameter.SeqNum=0;
    CParameter.UniqueSeqNum=0;
    CParameter.PrevIterk=0;
    CParameter.FinalCenterID=NULL;
    CParameter.FinalCenterSeq=NULL;
    CParameter.FinalSD=NULL;
    CParameter.FinalSize=NULL;
    CParameter.FinalMembership=NULL;
    CParameter.r=NULL;
}

bool ChompFile(const char* fname, int block, CROPParameters &CParameter)
{
    fstream in,out;
    string ResultName;
    string ResultNamePrefix(fname);
    ResultNamePrefix+=".BlockSeq";
    string ResultNameNum="1";
    string Line;
    char Buffer[ReadSize];
    char temp[ReadSize];
    in.open(fname);
    if(in.fail())
    {
        cout<<"Read Fasta File Error"<<endl;
        return 0;
    }
    while(!in.eof())
    {
        in.getline(temp, ReadSize);
        if(temp[0]=='>')
            CParameter.UniqueSeqNum++;
    }
    if(CParameter.UniqueSeqNum==0)
    {
        cout<<"Not A Valid Fasta File!"<<endl;
        return 0;
    }    
    in.clear();   
    in.seekg(0,ios::beg);
    cout<<CParameter.UniqueSeqNum<<endl;
    ResultName=ResultNamePrefix+ResultNameNum; 
    out.open(ResultName.c_str(),ios_base::out);
    if(out.fail())
    {
        cout<<"Write Fasta File Error"<<endl;
    }
    in.getline(temp,ReadSize);
    Line.assign(temp);
    if(temp[0]!='>')
    {
        cout<<"Bad Header"<<endl;
        return 0;
    }
    out<<Line<<endl;
    in.getline(temp,ReadSize);
    Line.assign(temp);
    int numblock=1;
    for(int i=0; i<CParameter.UniqueSeqNum; i++)
    {
        while(temp[0]!='>')
        {
            out<<Line<<endl;
            in.getline(temp,ReadSize);
            if(temp[0] == '\n' || temp[0] == '\0') break;
            Line.assign(temp);
        }
        if(floor((double)MAX(i-1,0)/CParameter.UniqueSeqNum*block)!=floor((double)i/CParameter.UniqueSeqNum*block))
        { 
            numblock++;
            out.close();
            memset(Buffer, 0, sizeof(char)*ReadSize);
            itoa(numblock, Buffer, 10);
            ResultNameNum.assign(Buffer);
            ResultName=ResultNamePrefix+ResultNameNum;
            out.open(ResultName.c_str(),ios_base::out);
            if(out.fail())
            {
                cout<<"Write Fasta File Error"<<endl;
            }
        }
        if(i<CParameter.UniqueSeqNum-1)
        {
            out<<Line<<endl;
            in.getline(temp,ReadSize);
            Line.assign(temp);
        }
        else
            out.close();
    }
    in.close();
    return 1;
}

bool ChompInterFile(const char* fname, int block, CROPParameters &CParameter)
{
    int tempstat=0;
    fstream in,out;
    string ResultName;
    string ResultNamePrefix(fname);
    ResultNamePrefix+=".BlockSeq";
    string ResultNameNum="1";
    string Line;
    int tempint;
    double tempdouble;
    char Buffer[ReadSize];
    char temp[ReadSize];
    string TempCenters(fname);
    TempCenters+=".TempCenters";
    in.open(TempCenters.c_str());
    if(in.fail())
    {
        cout<<"Read Inter Result Fasta File Error"<<endl;
        return 0;
    }
    while(!in.eof())
    {
        in.getline(temp, ReadSize);
        if(temp[0]=='>')
            tempstat++;
    }
    if(tempstat==0)
    {
        cout<<"Not A Valid Fasta File!"<<endl;
        return 0;
    }    
    in.clear();   
    in.seekg(0,ios::beg);
    cout<<"Number of Seqs in File Before Chomp = "<<tempstat<<endl;
    ResultName=ResultNamePrefix+ResultNameNum; 
    out.open(ResultName.c_str(),ios_base::out);
    int numblock=1;
    if(out.fail())
    {
        cout<<"Write Fasta File Error"<<endl;
    }
    for(int i=0;i<CParameter.Finalk;i++)
    {
        if(floor((double)MAX(i-1,0)/CParameter.Finalk*block)!=floor((double)i/CParameter.Finalk*block))
        {
            numblock++;
            out.close();
            itoa(numblock, Buffer, 10);
            ResultNameNum.assign(Buffer);
            ResultName=ResultNamePrefix+ResultNameNum;
            out.open(ResultName.c_str(),ios_base::out);
            if(out.fail())
            {
                cout<<"Write Fasta File Error"<<endl;
            }
        }
        //*******Center Seq ID**********
        in.getline(temp,ReadSize);
        Line.assign(temp);
        out<<Line<<endl;
        //*******Center Seq************
        in.getline(temp,ReadSize);
        Line.assign(temp);
        out<<Line<<endl;
        //********Abudance and Sigma********
        in>>tempint;
        out<<tempint<<endl;
        in>>tempdouble;
        out<<tempdouble<<endl;
        in.getline(temp,ReadSize);
        //***********Eigen Seqs**************
        for(int j=0;j<MIN(tempint,20);j++)
        {
            in.getline(temp,ReadSize);
            Line.assign(temp);
            out<<Line<<endl;
        }
        if(i==CParameter.Finalk-1)
            out.close();
    }
    in.close();
    return 1;
}

int LoadingList(const char* fname, CROPParameters &CParameter)
{
    int N;
    fstream SeqList;
    string SeqListName(fname);
    char dummy[MAX_LINE];
    SeqListName+=".list";
    cout<<"Loading "<<SeqListName<<endl;
    SeqList.open(SeqListName.c_str());
    if(SeqList.fail())
    {
        cout<<"Cannot open name file"<<endl;
        return 0;
    }
    SeqList>>N;
    SeqList.getline(dummy, MAX_LINE);
    for(int i=0;i<N;i++)
    {
        string temp1, temp2, temp;
        int len = 0;
        SeqList>>temp1;
        SeqList.getline(dummy, MAX_LINE);//read .Name
        //len = strlen(dummy);
        //len += (dummy[len-1] == '\n')? -1:0;
        //temp1.assign(dummy,len);
        SeqList>>temp2;
        SeqList.getline(dummy, MAX_LINE);//read .Members
        len = strlen(dummy);        
        CParameter.List[temp1]=temp2;  //build map
    }
    SeqList.close();
    return N;
}

/*
void MergeResult(int a)
{
    bool flag;
    if(a==0)
    {
        FinalCenterID=new string[SeqNum];
        FinalCenterSeq=new string[SeqNum];
        FinalSize=new int[SeqNum];
        FinalSD=new double[SeqNum];
    }
    for(int i=0;i<BlockResult[a].ClusterNumber;i++)
    {
            flag=true;
            for(int j=0;j<Finalk;j++)
            {
                int matrix = EDNAFULL;
           	    alignment (*alignTool) (const char *, const char *, alignment, int, int) = NULL, align = {0, 0, 0};
                alignTool = NW_alignmentEx;
			    align = alignTool(FinalCenterSeq[j].c_str(), BlockResult[a].ClusterCenterSeq[i], align, matrix, 1);
			    double d=(1.0-(double)align.matches/(double)SHORT_ALIGN(align))*100.0;
		        if(d>15)
		            continue;
                else
                {
                    int m=BlockResult[a].ClusterSize[i];
                    int n=FinalSize[j];
                    double SD1=BlockResult[a].ClusterStandardDeviation[i];
                    double SD2=FinalSD[j];
                    double Y1=(SD1*SD1*m+(SD1+d/2)*(SD1+d/2)*n)/(m+n);
                    double Y2=((SD1+d/2)*(SD1+d/2)*m+SD2*SD2*n)/(m+n);
                    double Y3=(SD1*SD1*m+SD2*SD2*n)/(m+n);
                    double MY=MIN(Y1,Y2);
                    if((Y1<(2+(int)ctrl*3)*(2+(int)ctrl*3))||(Y2<(2+(int)ctrl*3)*(2+(int)ctrl*3)))
                    //if((sqrt((double)(m*m+n*n)/(m+n)/(m+n)*(pow(SD1,4)+pow(SD2,4)))<(2+(int)ctrl*3)*(2+(int)ctrl*3)))
                    {
                        if(MY==Y1)
                            FinalCenterSeq[j].assign(BlockResult[a].ClusterCenterSeq[i]);
                        FinalSize[j]+=m;
                        FinalSD[j]=sqrt(Y3);
                        flag=false;
                        cout<<"Merge ["<<a<<"]["<<i<<"] To ["<<j<<"]"<<endl;
                        break;
                    }
                    else
                        continue;
                }
            }
            if(flag==true)
            {
                FinalCenterID[Finalk].assign(BlockResult[a].ClusterCenterID[i]);
                FinalCenterSeq[Finalk].assign(BlockResult[a].ClusterCenterSeq[i]);
                FinalSize[Finalk]=BlockResult[a].ClusterSize[i];
                FinalSD[Finalk]=BlockResult[a].ClusterStandardDeviation[i];
                Finalk++;
                cout<<"Add ["<<a<<"]["<<i<<"]"<<endl;
            }
    }
//**********final refine***************
    if(a==Block-1){
    for(int i=Finalk-1;i>-1;i--)
    {
            for(int j=0;j<i;j++)
            {
                int matrix = EDNAFULL;
           	    alignment (*alignTool) (const char *, const char *, alignment, int, int) = NULL, align = {0, 0, 0};
                alignTool = NW_alignmentEx;
			    align = alignTool(FinalCenterSeq[j].c_str(), FinalCenterSeq[i].c_str(), align, matrix, 1);
			    double d=(1.0-(double)align.matches/(double)LONG_ALIGN(align))*100.0;
		        if(d>15)
		            continue;
                else
                {
                    int m=FinalSize[i];
                    int n=FinalSize[j];
                    double SD1=FinalSD[i];
                    double SD2=FinalSD[j];
                    double Y1=(SD1*SD1*m+(SD1+d/2)*(SD1+d/2)*n)/(m+n);
                    double Y2=((SD1+d/2)*(SD1+d/2)*m+SD2*SD2*n)/(m+n);
                    double MY=MIN(Y1,Y2);
                    if((Y1<(2+(int)ctrl*3)*(2+(int)ctrl*3))||(Y2<(2+(int)ctrl*3)*(2+(int)ctrl*3)))
                    //if((sqrt((double)(m*m+n*n)/(m+n)/(m+n)*(pow(SD1,4)+pow(SD2,4)))<(2+(int)ctrl*3)*(2+(int)ctrl*3)))
                    {
                        if(MY==Y1)
                            FinalCenterSeq[j]=FinalCenterSeq[i];
                        FinalSize[j]+=m;
                        FinalSD[j]=sqrt(MY);
                        cout<<"Merge ["<<i<<"] To ["<<j<<"]"<<endl;
                        Finalk--;
                        break;
                    }
                    else
                        continue;
                }
            }
    }
    }
    return;
}
*/

void MergeUsingBC(const char* fname, CROPParameters &CParameter)
{
    cout<<"MergingUsingBC"<<endl;
    string TempCenters(fname);
    TempCenters+=".TempCenters";
    BayesianClusteringResult MergeCenterResult;
    Parameters Parameter;
    BayesianClustering(TempCenters.c_str(),CParameter.Finalk*10, CParameter.Lower, CParameter.Upper, MergeCenterResult, 20, Parameter);
    //cout<<"Bayesian Clustering Completed"<<endl;
    CParameter.Finalk=MergeCenterResult.ClusterNumber;
    CParameter.FinalCenterID=new string[CParameter.Finalk];
    CParameter.FinalCenterSeq=new string[CParameter.Finalk];
    CParameter.FinalSize=new int[CParameter.Finalk];
    CParameter.FinalSD=new double[CParameter.Finalk];
    CParameter.FinalMembership=new string*[CParameter.Finalk];
    string** TempFinalMembership;
    TempFinalMembership=new string*[CParameter.Finalk];
    memset(CParameter.FinalSize,0,sizeof(int)*CParameter.Finalk);
    for(int i=0;i<CParameter.Finalk;i++)
    {
        CParameter.FinalCenterID[i].assign(MergeCenterResult.ClusterCenterID[i]);
        CParameter.FinalCenterSeq[i].assign(MergeCenterResult.ClusterCenterSeq[i]);
        CParameter.FinalSD[i]=MergeCenterResult.ClusterStandardDeviation[i];
        TempFinalMembership[i]=new string[CParameter.SeqNum];
        for(int p=0;p<MergeCenterResult.ClusterIterationSize[i];p++)
        {
            for(int j=0;j<CParameter.Block[CParameter.Maxiter];j++)
            {
                for(int l=0;l<CParameter.BlockResult[CParameter.Maxiter][j].ClusterNumber;l++)
                {
                    if(strcmp(MergeCenterResult.ClusterMember[i][p].c_str(),CParameter.BlockResult[CParameter.Maxiter][j].ClusterCenterID[l])==0)
                    {                                                                                                                                                                
                            for(int q=0;q<CParameter.BlockResult[CParameter.Maxiter][j].ClusterIterationSize[l];q++)
                            {
                                if(CParameter.BlockResult[CParameter.Maxiter][j].ClusterMemberFlag[l][q]==true)
                                {
                                    cout<<"Redundant Sequence Detected: "<<CParameter.BlockResult[CParameter.Maxiter][j].ClusterMember[l][q]<<endl;
                                }
                                else
                                {
                                    CParameter.BlockResult[CParameter.Maxiter][j].ClusterMemberFlag[l][q]=true;
                                    TempFinalMembership[i][CParameter.FinalSize[i]]=CParameter.BlockResult[CParameter.Maxiter][j].ClusterMember[l][q];
                                    CParameter.FinalSize[i]++;                                    
                                }
                            }//q
                    }//if
                }//l
            }//j
        }//p
    }//i
    int totalsize1,totalsize2;
    for(int i=0;i<CParameter.Maxiter+1;i++)
    {
        totalsize1=0;
        totalsize2=0;
        for(int j=0;j<CParameter.Block[i];j++)
            totalsize1+=CParameter.BlockResult[i][j].ClusterNumber;
        //cout<<"Total Number of Clusters "<<totalsize1<<endl;
        totalsize1=0;
        for(int j=0;j<CParameter.Block[i];j++)
            for(int l=0;l<CParameter.BlockResult[i][j].ClusterNumber;l++)
            {
                totalsize1+=CParameter.BlockResult[i][j].ClusterIterationSize[l];
                totalsize2+=CParameter.BlockResult[i][j].ClusterSize[l];
            }
       //cout<<"Total Seq Nums in Clusters "<<totalsize2<<" ,While Iteration Seq Nums Is "<<totalsize1<<endl;  
    }
    totalsize1=0;
    for(int i=0;i<CParameter.Finalk;i++)
        totalsize1+=CParameter.FinalSize[i];
    //cout<<"Currently Assigned Seq Num is "<<totalsize1<<endl;
    /*
    for(int i=0;i<Block[Maxiter];i++)
        for(int j=0;j<BlockResult[Maxiter][i].ClusterNumber;j++)
            for(int l=0;l<BlockResult[Maxiter][i].ClusterIterationSize[j];l++)
                if(BlockResult[Maxiter][i].ClusterMemberFlag[j][l]==false)
                    cout<<"Missing Sequence Detected: "<<BlockResult[Maxiter][i].ClusterMember[j][l]<<endl;
    */
    int TempFinalSize[CParameter.Finalk];
    while(CParameter.Maxiter>0)
    {
        CParameter.Maxiter--;
        memcpy(TempFinalSize, CParameter.FinalSize, sizeof(int)*CParameter.Finalk);
        for(int i=0;i<CParameter.Finalk;i++)
        {
            for(int p=0;p<TempFinalSize[i];p++)
            {
                for(int j=0;j<CParameter.Block[CParameter.Maxiter];j++)
                {
                    for(int l=0;l<CParameter.BlockResult[CParameter.Maxiter][j].ClusterNumber;l++)
                    {
                       if(strcmp(TempFinalMembership[i][p].c_str(),CParameter.BlockResult[CParameter.Maxiter][j].ClusterCenterID[l])==0)
                       {
                           for(int q=0;q<CParameter.BlockResult[CParameter.Maxiter][j].ClusterIterationSize[l];q++)
                           {
                               if(strcmp(CParameter.BlockResult[CParameter.Maxiter][j].ClusterMember[l][q].c_str(),CParameter.BlockResult[CParameter.Maxiter][j].ClusterCenterID[l])!=0)
                               {
                                    if(CParameter.BlockResult[CParameter.Maxiter][j].ClusterMemberFlag[l][q]==true)
                                    {
                                        cout<<"Redundant Sequence Detected: "<<CParameter.BlockResult[CParameter.Maxiter][j].ClusterMember[l][q]<<endl;
                                    }
                                    else
                                    {
                                        CParameter.BlockResult[CParameter.Maxiter][j].ClusterMemberFlag[l][q]=true;
                                        TempFinalMembership[i][CParameter.FinalSize[i]]=CParameter.BlockResult[CParameter.Maxiter][j].ClusterMember[l][q];
                                        CParameter.FinalSize[i]++; 
                                    }
                               }//if
                           }//q
                       }//if
                   }//l
                }//j
            }//p
        }//i
        totalsize1=0;
        for(int i=0;i<CParameter.Finalk;i++)
         totalsize1+=CParameter.FinalSize[i];
        //cout<<"Currently Assigned Seq Num is "<<totalsize1<<endl;
   }//while
   memcpy(TempFinalSize, CParameter.FinalSize, sizeof(int)*CParameter.Finalk);
   for(int i=0; i<CParameter.Finalk; i++)
   {
       for(int j=0; j<TempFinalSize[i];j++)
       {
           char* pch;
           char* temp;
           int len = strlen(CParameter.List[TempFinalMembership[i][j]].c_str());
           temp=new char[len+1];
           memset(temp, 0x00, sizeof(char)*(len+1));
           strcpy(temp, CParameter.List[TempFinalMembership[i][j]].c_str());
           pch = strtok (temp, ",");
           pch = strtok (NULL, ",");
           while (pch != NULL)
           {
                 TempFinalMembership[i][CParameter.FinalSize[i]].assign(pch);
                 CParameter.FinalSize[i]++;
                 pch = strtok (NULL, ",");
           }
           delete[] pch;
           delete[] temp;
       }
   }
   //cout<<"reading list file finished"<<endl;
   totalsize1=0;
   for(int i=0;i<CParameter.Finalk;i++)
       totalsize1+=CParameter.FinalSize[i];
   //cout<<"Currently Assigned Seq Num is "<<totalsize1<<endl; 
   for(int i=0;i<CParameter.Finalk;i++)
       if(CParameter.FinalSize[i]!=MergeCenterResult.ClusterSize[i])
           cout<<CParameter.FinalSize[i]<<" "<<MergeCenterResult.ClusterSize[i]<<" Unexpected Merging Result!"<<endl;       
   for(int i=0;i<CParameter.Finalk;i++)
   {
       CParameter.FinalMembership[i]=new string[CParameter.FinalSize[i]];
       for(int j=0;j<CParameter.FinalSize[i];j++)
           CParameter.FinalMembership[i][j]=TempFinalMembership[i][j];
   }
   
   //*********Free Memory*****************
   for(int i=0; i<CParameter.Finalk; i++)
   {
        delete[] MergeCenterResult.ClusterMember[i];
        delete[] MergeCenterResult.ClusterCenterID[i];
        delete[] MergeCenterResult.ClusterCenterSeq[i];
        delete[] TempFinalMembership[i];
   }
   delete[] MergeCenterResult.ClusterMember;
   delete[] MergeCenterResult.ClusterCenterID;
   delete[] MergeCenterResult.ClusterCenterSeq;
   delete MergeCenterResult.ClusterStandardDeviation;
   delete MergeCenterResult.ClusterSize;
   delete[] TempFinalMembership;
   //**************************************
                                              
   return;                      
}

void InterResultRead(const char* fname, int iter, CROPParameters &CParameter)
{
    string TempCenters(fname);
    TempCenters+=".TempCenters";
    fstream BlockCenters;
    CParameter.Finalk=0;
    BlockCenters.open(TempCenters.c_str(),ios_base::out);
    if(BlockCenters.fail())
    {
        cout<<"Fail to Open Temp File"<<endl;
        return;
    }
    int totalcluster=0;
    for(int i=0;i<CParameter.Block[iter];i++)
        for(int j=0;j<CParameter.BlockResult[iter][i].ClusterNumber;j++)
            totalcluster++;
    int* Shuffle;
    Shuffle=new int[totalcluster];
    int Count=0;
    for(int i=0;i<CParameter.Block[iter];i++)
        for(int j=0;j<CParameter.BlockResult[iter][i].ClusterNumber;j++)
        {
            Shuffle[Count]=i;
            Count++;
        }
    CParameter.r = gsl_rng_alloc (gsl_rng_taus);
    gsl_rng_set(CParameter.r,(unsigned int)time((time_t *)NULL));
    gsl_ran_shuffle(CParameter.r,Shuffle,totalcluster,sizeof(int));
    int TempCount[CParameter.Block[iter]];
    memset(TempCount,0,sizeof(int)*CParameter.Block[iter]);
    for(int i=0;i<totalcluster;i++)
    {
        BlockCenters<<">"<<CParameter.BlockResult[iter][Shuffle[i]].ClusterCenterID[TempCount[Shuffle[i]]]<<"\n";
        BlockCenters<<CParameter.BlockResult[iter][Shuffle[i]].ClusterCenterSeq[TempCount[Shuffle[i]]]<<"\n";
        BlockCenters<<CParameter.BlockResult[iter][Shuffle[i]].ClusterSize[TempCount[Shuffle[i]]]<<"\n";
        BlockCenters<<CParameter.BlockResult[iter][Shuffle[i]].ClusterStandardDeviation[TempCount[Shuffle[i]]]<<"\n";
        for(int l=0;l<MIN(20,CParameter.BlockResult[iter][Shuffle[i]].ClusterSize[TempCount[Shuffle[i]]]);l++)
        {
            BlockCenters<<CParameter.BlockResult[iter][Shuffle[i]].ClusterMemberSeq[TempCount[Shuffle[i]]][l];
            if((i!=CParameter.Block[iter]-1)||(TempCount[Shuffle[i]]!=CParameter.BlockResult[iter][Shuffle[i]].ClusterNumber-1)||(l!=MIN(20,CParameter.BlockResult[iter][Shuffle[i]].ClusterSize[TempCount[Shuffle[i]]])-1))
                BlockCenters<<"\n";
        }//l
        CParameter.Finalk++;        
        TempCount[Shuffle[i]]++;
    }//i
    BlockCenters.close();
    delete[] Shuffle;
    gsl_rng_free(CParameter.r);
    cout<<"Number of Clusters in Previous Round = "<<CParameter.Finalk<<endl;
    return;
}

int CROP(string fname, string oname, int TempBlockSize, int TempStep, int TempCtrl, float Lower, float Upper, int TempMaxiter, int TempMaxiterBlockSize, CROPParameters &CParameter)
{
    time_t start,stop;
    fstream finalout;
    fstream logfile;
    string logfilename;
    logfilename=oname+".log";
    logfile.open(logfilename.c_str(), ios_base::out);
    if(logfile.fail())
    { 
        cout<<"cannot create log file"<<endl;
        return 0;
    }
	start=time(NULL);
	InitCROPParameters(CParameter);
	CParameter.Block[0]=TempBlockSize;
	CParameter.step=TempStep;
    switch(TempCtrl)
    {
        case 0:
             CParameter.Lower=1.5; 
             CParameter.Upper=2.5;
             break;
        case 1:
             CParameter.Lower=1;
             CParameter.Upper=1.5;
             break;
        case 2:
             CParameter.Lower=Lower;
             CParameter.Upper=Upper;
             break;
        default:
             cout<<"no such threhold options, -g will be used"<<endl;
             logfile<<"-g is used\n";
             CParameter.Lower=1.5; 
             CParameter.Upper=2.5;
             break;
    }
	CParameter.Maxiter=TempMaxiter;
	CParameter.MaxiterBlockSize=TempMaxiterBlockSize;
	CParameter.BlockResult=new BayesianClusteringResult*[CParameter.Maxiter];
	CParameter.BlockResult[0]=new BayesianClusteringResult[CParameter.Block[0]];
	CParameter.SeqNum=LoadingList(fname.c_str(),CParameter);
	logfile<<"-i "<<fname<<" -b "<<TempBlockSize<<" -e "<<TempStep<<" -m "<<TempMaxiter<<" -z "<<TempMaxiterBlockSize<<" Threshold: "<<CParameter.Lower<<" "<<CParameter.Upper<<endl;
	logfile<<CParameter.SeqNum<<endl;
    if(!ChompFile(fname.c_str(), CParameter.Block[0], CParameter))
    {
        cout<<"chop file fail"<<endl;
        logfile<<"chop file fail"<<endl;
        return 0;
    }
    #pragma omp parallel for
    for(int i=0;i<CParameter.Block[0];i++)
    {
        char Buf[ReadSize];
        memset(Buf, 0, sizeof(char)*ReadSize);
        Parameters Parameter;
        itoa(i+1,Buf,10);
        string BlockFnameNum(Buf);
        string BlockFnameString=fname;
        BlockFnameString+=".BlockSeq"+BlockFnameNum;
        #pragma omp critical
        {
            cout<<"Round: 0 Block: "<<i+1<<endl;
        }//omp critical
        BayesianClustering(BlockFnameString.c_str(), CParameter.step, CParameter.Lower, CParameter.Upper, CParameter.BlockResult[0][i],0,Parameter);
        #pragma omp critical
        {
            logfile<<"Round: 0 Block: "<<i+1<<": "<<CParameter.BlockResult[0][i].ClusterNumber<<endl;
        }
    } 
    for(int i=1;i<CParameter.Maxiter;i++)
    {
        InterResultRead(fname.c_str(), i-1, CParameter);
        if((CParameter.Finalk>=(double)0.9*CParameter.PrevIterk)&&(i>1))
        {
            CParameter.Maxiter=i-1;                 
            break;
        }//if
        if(CParameter.Finalk>CParameter.MaxiterBlockSize)
        {
            //2010-2-7 if merging is not efficient, try to increase the block size
            //********************************************************************
            if(CParameter.Finalk>=(double)0.8*CParameter.PrevIterk)
            {
                CParameter.PrevIterk=CParameter.Finalk;
                CParameter.Block[i]=MAX(CParameter.Finalk/MIN(1.5*CParameter.MaxiterBlockSize, LIMITB)+1,omp_get_num_procs());
            }
            else
            {
                CParameter.PrevIterk=CParameter.Finalk;
                CParameter.Block[i]=MAX(CParameter.Finalk/CParameter.MaxiterBlockSize+1,omp_get_num_procs());
            }
            //********************************************************************
            //2010-2-7 done
            if(!ChompInterFile(fname.c_str(), CParameter.Block[i], CParameter))
            {
                cout<<"Chomp Fail!"<<endl;
                logfile<<"Chomp Fail!"<<endl;
                return 0;
            }
            CParameter.BlockResult[i]=new BayesianClusteringResult[CParameter.Block[i]];
            #pragma omp parallel for
            for(int j=0;j<CParameter.Block[i];j++)
            {
                char Buf[ReadSize];
                Parameters Parameter;
                itoa(j+1,Buf,10);
                string BlockFnameNum(Buf);
                string BlockFnameString=fname;
                BlockFnameString+=".BlockSeq"+BlockFnameNum;
                #pragma omp critical
                {
                    cout<<"Round: "<<i<<" Block: "<<j+1<<endl;
                }//omp critical
                BayesianClustering(BlockFnameString.c_str(), CParameter.step, CParameter.Lower, CParameter.Upper, CParameter.BlockResult[i][j],i,Parameter);
                #pragma omp critical
                {
                    logfile<<"Round: "<<i<<" Block: "<<j+1<<": "<<CParameter.BlockResult[i][j].ClusterNumber<<endl;
                }
            }//j
            if(i==CParameter.Maxiter-1)
            {
                InterResultRead(fname.c_str(), i, CParameter);
                CParameter.Maxiter=i;
            }
        }//if
        else
        {
            CParameter.Maxiter=i-1;
            break;
        }//else
    }//i
    //finalout<<endl;
    //finalout<<"Merging Using BC"<<endl;
    logfile<<"Merging Using BC"<<endl;
    MergeUsingBC(fname.c_str(), CParameter);
    logfile<<"Writing Output"<<endl;
    logfile.close();
    /*
    for(int i=0;i<Block;i++)
        MergeResult(i);
    */       
    //***************output result***************************************
    string ofilename;
    ofilename=oname+".cluster";
    finalout.open(ofilename.c_str(), ios_base::out);
    if(finalout.fail())
    { 
        cout<<"cannot open output file 1"<<endl;
        return 0;
    }
    finalout<<CParameter.Finalk<<endl;
    for(int i=0;i<CParameter.Finalk;i++)
    {
        finalout<<CParameter.FinalCenterID[i]<<"\t"<<CParameter.FinalSize[i]<<"\t"<<CParameter.FinalSD[i]<<"\n";
    }
    finalout.close();
    ofilename=oname+".cluster.fasta";
    finalout.open(ofilename.c_str(), ios_base::out);
    if(finalout.fail())
    { 
        cout<<"cannot open output file 2"<<endl;
        return 0;
    }
    for(int i=0;i<CParameter.Finalk;i++)
    {
        finalout<<">"<<CParameter.FinalCenterID[i]<<"\n";
        finalout<<CParameter.FinalCenterSeq[i]<<"\n";
    }
    finalout.close();
    ofilename=oname+".cluster.list";
    finalout.open(ofilename.c_str(), ios_base::out);
    if(finalout.fail())
    { 
        cout<<"cannot open output file 3"<<endl;
        return 0;
    }
    for(int i=0;i<CParameter.Finalk;i++)
    {
        finalout<<CParameter.FinalCenterID[i]<<"\t";
        for(int j=0;j<CParameter.FinalSize[i];j++)
        {
            finalout<<CParameter.FinalMembership[i][j];
            if(j==CParameter.FinalSize[i]-1)
                finalout<<"\n";
            else
                finalout<<",";
        }
    }
    ofilename.clear();
    finalout.close();
    //cout<<"output done"<<endl;
    //******************************************************************
    //**********Free Memory*********************************************
    DeleteTmpFiles(fname, CParameter);
    for(int i=0;i<CParameter.Maxiter+1;i++)
    {
        for(int j=0;j<CParameter.Block[i];j++)
        {
            for(int l=0; l<CParameter.BlockResult[i][j].ClusterNumber; l++)
            {
                delete[] CParameter.BlockResult[i][j].ClusterMember[l];
                delete[] CParameter.BlockResult[i][j].ClusterMemberFlag[l];
                delete[] CParameter.BlockResult[i][j].ClusterCenterID[l];
                delete[] CParameter.BlockResult[i][j].ClusterCenterSeq[l];
                delete[] CParameter.BlockResult[i][j].ClusterMemberSeq[l];
            }//l
            delete[] CParameter.BlockResult[i][j].ClusterMemberSeq;
            delete[] CParameter.BlockResult[i][j].ClusterMember;
            delete[] CParameter.BlockResult[i][j].ClusterMemberFlag;
            delete[] CParameter.BlockResult[i][j].ClusterCenterID;
            delete[] CParameter.BlockResult[i][j].ClusterCenterSeq;
            delete[] CParameter.BlockResult[i][j].ClusterStandardDeviation;
            delete[] CParameter.BlockResult[i][j].ClusterSize;
            delete[] CParameter.BlockResult[i][j].ClusterIterationSize;
        }//j
        delete[] CParameter.BlockResult[i];
    }
    delete[] CParameter.BlockResult;
    delete[] CParameter.FinalCenterID;
    delete[] CParameter.FinalCenterSeq;
    delete[] CParameter.FinalSD;
    delete[] CParameter.FinalSize;
    for(int i=0;i<CParameter.Finalk;i++)
        delete[] CParameter.FinalMembership[i];
    delete[] CParameter.FinalMembership;
    CParameter.List.empty();
    //**********************************************
    stop=time(NULL);
    cout<<"Time Cost Is "<<stop-start<<" Sec"<<endl;
    return 1;    
}

void DeleteTmpFiles(string fname, CROPParameters &CParameter)
{
    char Buffer[10];
    for(int i=1;i<=CParameter.Block[0];i++)
    {
        string tempNum(itoa(i, Buffer, 10));
        string tempFile=fname+".BlockSeq"+tempNum;
        if(remove(tempFile.c_str())!=0)
            cout<<"Error deleting "<<tempFile<<endl;
    }
    string tempFile=fname+".TempCenters";
    if(remove(tempFile.c_str())!=0)
        cout<<"Error deleting "<<tempFile<<endl;
    return;
} 
